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H , Abstract 

We examine a (3 + l)-dimensional model of staggered lattice fermions 
with a four-fermion interaction and Z(2) chiral symmetry using the 
Hamiltonian formulation. This model cannot be simulated with stan- 
dard fermion algorithms because those suffer from a very severe sign 
problem. We use a new fermion simulation technique — the meron- 
cluster algorithm — which solves the sign problem and leads to high- 
precision numerical data. We investigate the finite temperature chiral 
phase transition and verify that it is in the universality class of the 3-d 
Ising model using finite-size scaling. 
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1 Introduction 

The numerical simulation of lattice fermions is a notoriously difficult problem which 
is the major stumbling block in solving QCD and other fermionic field theories. 
The standard method is to integrate out the fermions and to simulate the resulting 
bosonic problem with a non-local action. In several cases of physical interest — for 
example, for QCD with an odd number of fiavors or with non-zero chemical potential 
— the bosonic Boltzmann factor may become negative or even complex and thus 
cannot be interpreted as a probability. When the sign or the complex phase of 
the Boltzmann factor is included in measured observables, the numerical simulation 
suffers from severe cancellations resulting in a sign problem. The standard fermion 
algorithms are incapable of exploring such models. As a consequence, QCD is usually 
simulated with an even number of flavors and at zero chemical potential. Even in the 
absence of a sign problem, the simulation of fermions is difficult. For example, lattice 
QCD simulations suffer from critical slowing down when one approaches the chiral 
limit in which the quarks become massless. In particular, this makes it difficult to 
identify the universality class of the flnite temperature QCD chiral phase transition. 

Even in simpler models with four-fermion interactions the identiflcation of the 
flnite temperature critical behavior is a non-trivial issue ^. A model with N fermion 
flavors shows mean-fleld behavior in the N = oo limit. On the other hand, at flnite 
N one flnds the non-trivial critical behavior that one expects based on dimensional 
reduction and standard universality arguments. For example, in it has been 
verifled that the chiral phase transition in a (2 -|- l)-d four-fermion model with 
A^ = 4 and Z(2) chiral symmetry is in the universality class of the 2-d Ising model. 
Due to the fermion sign problem, standard fermion simulation methods often do not 
work in models with a too small number of flavors. 

In this paper we apply a recently developed technique for solving the sign 
problem to a (3 -|- l)-d model of staggered fermions using the Hamiltonian formula- 
tion. The model has N = 2 flavors and a Z(2) chiral symmetry that is spontaneously 
broken at low temperatures. The fermion determinant can be negative in this model. 
Hence, due to the sign problem standard fermion algorithms fail in this case. Our 
algorithm is the only numerical method available to simulate this model. In this 
method we do not integrate out the fermions but describe them in a Fock state basis. 
The resulting bosonic model of fermion occupation numbers interacts locally, but 
has a non-local fermion permutation sign resulting from the Pauli exclusion princi- 
ple. Standard numerical methods would suffer from severe cancellations of positive 
and negative contributions to the partition function. Like other cluster methods, our 
algorithm decomposes a conflguration of fermion occupation numbers into clusters 
which can be flipped independently. Under a cluster flip an occupied site becomes 
empty and vice versa. The main idea of the meron-cluster algorithm is to construct 
the clusters such that they affect the fermion sign independent of each other when 
they are flipped. In addition, it must always be possible to flip the clusters into a ref- 



erence configuration with a positive sign. A cluster whose flip changes the fermion 
sign is referred to as a meron because it can be viewed as a half-instanton. If a 
conflguration contains a meron-cluster, its contribution to the partition function is 
canceled by the contribution of the conflguration that one obtains when the meron- 
cluster is flipped. The observables that we consider get non-zero contributions from 
the zero- and two-meron sectors only. Our algorithm ensures that conflgurations 
with more than two merons are never generated, which leads to an exponential gain 
in statistics and to a complete solution of the sign problem. 

Like other cluster algorithms the meron algorithm substantially reduces critical 
slowing down. This allows us to work directly in the chiral limit. As a result, we 
can study the nature of the chiral phase transition in great detail. The Z(2) chiral 
symmetry is spontaneously broken at low temperatures and gets restored in the 
high-temperature phase. As expected, the system close to the flnite temperature 
critical point is in the universality class of the 3-d Ising model. We verify this in a 
high-precision flnite-size scaling investigation of the chiral susceptibility. 

The paper is organized as follows. In section 2 we introduce the staggered fermion 
Hamiltonian, derive the path integral representation of the partition function, and 
discuss the fermion sign as well as relevant observables. Section 3 contains the de- 
scription of the meron-cluster algorithm and the corresponding improved estimators. 
In section 4, we present the results of numerical simulations and extract the critical 
behavior from the flnite-size scaling of the chiral susceptibility. Finally, section 5 
contains our conclusions. 



2 The Staggered Fermion Model 

Let us consider staggered fermions hopping on a 3-d cubic spatial lattice with V = L^ 
sites X {L even) and with periodic or antiperiodic spatial boundary conditions. We 
start in the Hamiltonian formulation and then derive a path integral on a (3 + l)-d 
Euclidean space-time lattice. The fermions are described by creation and annihila- 
tion operators \E'^ and "^x with standard anticommutation relations 

{M/+, M/ + } = {M/., ^y} = 0, {Vi/+, ^y} = 6xy. (2.1) 

The staggered fermion Hamilton operator takes the form 

H = J2 h., + m^(-l)^i+^=^+^^vi/+vl/,, (2.2) 

x,i X 

that is a sum of nearest-neighbor couplings h^^i and a mass term m\E'\E'. In the 
following we work directly in the chiral limit, -m = 0, and only use \E'\E' as an 
observable. The term h^ i couples the fermion operators at the lattice sites x and 



X + i, where i is a unit-vector in the i-direction, and 

K, = Iv.A'^t^.^i + K^^ + Gi^t"^. - ^)(^:+.^.+i - 1)- (2.3) 

Here rj^^i = 1, rix2 = (~1)^^ and rj^^^ = (—1)^1+^2 g^j^g i\^q standard staggered fermion 
sign factors, and G is a four-fermion couphng constant. The system has a conserved 
fermion number 

N = J2^t^,, (2.4) 

X 

because [H, N] = 0. Besides the f/(l) fermion number symmetry, the model has a 
Z(2) chiral symmetry, which (up to a phase) simply shifts \1/^ and "^x by one lattice 
spacing in all three directions. This changes the sign of \I^\I^ but leaves the m = 
Hamiltonian invariant. There are also other Z(2) symmetries which correspond 
to discrete flavor transformations. For a detailed discussion of the symmetries of 
staggered fermions in the Hamiltonian formulation we refer to ^. 

To construct a path integral for the partition function, we decompose the Hamil- 
ton operator into six terms H = Hi + H2 + ■■■ + Hq, with 

Hi = ^ hx^i, Hi+3 = ^ hx^i- (2.5) 

x^{xi,X2,x^) x = (xi,X2,x^) 

XiCven Xiodd 

The individual contributions to a given Hi commute with each other, but two differ- 
ent Hi do not commute. Using the Suzuki- Trotter formula we express the fermionic 
partition function at inverse temperature (3 as 

Z. = Tr[exp(-/?i7)] = lim Tr[exp(-ei7i)exp(-ei72)-.exp(-ei/6)]^- (2.6) 

We have introduced 6M Euclidean time slices with e = (3/M being the lattice spacing 
in the Euclidean time direction. Following Jordan and Wigner we represent the 
fermion operators by Pauli matrices 

^+ = alal..atia+, ^, = alal..aUa^ , n, = ^+^, = -(af + 1), (2.7) 



with 



4 = I 

' 2 



.(^l ± ^o-f )> K, o-,^] = 2i6i„^eijk(Tf. (2. 



Here / labels the lattice point x. The Jordan- Wigner representation requires an 
ordering of the lattice points. For example, one can label the point x = {xi,X2,X3) 
(with Xi = 0, 1, ..., L — 1) by / = 1 + Xi + X2L + x^L^. It should be pointed out that 
the Jordan- Wigner representation works in any dimension. In one dimension the 
lattice points are, of course, naturally ordered, but even in higher dimensions the 
physics is completely independent of the arbitrary ordering. We now insert complete 
sets of fermion Fock states between the factors exp(— eifj). Each site is either empty 
or occupied, i.e. rix has eigenvalue or 1. In the Pauli matrix representation this 



corresponds to eigenstates |0) and |1) of af with o-f\0) = — 10) and crf\l) = |1). The 
transfer matrix is a product of factors 



eG 

exp(-eh^^i) = exp(— ) 



/ exp(-f ) 
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cosh 1 


S sinh f 





S sinh f 


cosh| 





V 





exp(-f ) ) 



(2.9) 



which is a 4 X 4 matrix in the Fock space basis |00), |01), |10) and |11) of two sites 
X and X + i. Here S = Vx,i<^f+i<^i+2- ■ ■'^m-i includes the local sign rj^^i as well as a 
non-local string of Pauli matrices running over consecutive labels between / and m, 
where / labels the lattice point x and m labels x + i. Note that S is diagonal in the 
occupation number basis. 

The partition function is now expressed as a path integral 

Zf = J2Sign[n]exp{-S[n]), (2.10) 

n 

over configurations of occupation numbers n{x,t) = 0, 1 on a (3 + l)-d space-time 
lattice of points {x,t). The Boltzmann factor 

exp{—S[n]) = J]^ exp{— s[n(a;, t),n(x + 1, t),n(x, t + l),n(x + 1, t + 1)]} 

xieven,t=6m 



X 



Y[ exp{— s[n(a;, t),n{x + 2, t),n{x, t + 1), n(x + 2, t + 1)]} 



x^(xi,X2,x^) 

X2even,t=6m+l 



X 



Y[ exp{—s[n{x, t),n{x + 3, t),n{x, t + 1), n(x + 3, t + 1)]} 



x^(xi,X2,x^) 

X3even,t=6m+2 



X 



Y[ exp{—s[n{x,t),n{x + i,t),n{x,t + l),n{x + 1, t + 1)]} 



x = (xi,X2,X3) 

xiodd,t=6m+3) 



X 



Y[ exp{— s[n(x, t),n{x + 2, t),n{x, t + l),n{x + 2, t + 1)]} 



x = (x-f^,X2,X3) 

X2odd,t=6m+A 



X 



Y[ exp{— s[n(a:, t),n{x + 3, t),n{x, t + l),n{x + 3, t + 1)]}, 



x = (x-f^,X2,X3) 

X3odd,t=6m+5 

(2.11) 
(with m = 0,l,...,M — 1) is a product of space-time plaquette contributions with 

exp(-s[0, 0, 0, 0]) = exp(-s[l, 1,1,1]) = exp(-^), 

exp(-s[0, 1, 0, 1]) = exp(-s[l, 0, 1,0]) = cosh |, 

exp(-s[0, 1,1,0]) = exp(-s[l, 0, 0, 1]) = sinh ^. (2.12) 



All the other Boltzmann weights are zero, which implies several constraints on al- 
lowed configurations. Note that here we have dropped the trivial overall factor 
exp(eG/4) that appeared in eq. (|2.9| ). 



The sign of a configuration, Sign[n], also is a product of space-time plaquette 
contributions sign[n(a:, t),n{x + i,t), n{x, t + 1), n(x + i^t + 1)] with 

sign[0, 0, 0, 0]) = sign[0, 1, 0, 1]) = sign[l, 0, 1, 0]) = sign[l, 1, 1, 1]) = 1, 
sign[0, 1, 1, 0]) = sign[l, 0, 0, 1]) = S. (2.13) 

It should be noted that S gets contributions from all lattice points with labels be- 
tween / and m. This seems to make an evaluation of the fermion sign rather tedious. 
Also, it is not a priori obvious that Sign[n] is independent of the arbitrarily chosen 
order of the lattice points. Fortunately, there is a simple way to compute Sign[ra], 
which is directly related to the Pauli exclusion principle and which is manifestly 
order-independent. In fact, Sign[n] has a topological meaning. The occupied lattice 
sites define fermion world-lines which are closed around the Euclidean time direc- 
tion. Of course, during their Euclidean time evolution fermions can interchange their 
positions, and the fermion world-lines define a permutation of particles. The Pauli 
exclusion principle dictates that the fermion sign is just the sign of that permuta- 
tion. If we work with antiperiodic spatial boundary conditions, Sign[n] receives an 
extra minus-sign for every fermion world-line that crosses a spatial boundary. Figure 
1 shows two configurations of fermion occupation numbers in (1 + 1) dimensions. 
The first configuration corresponds to two fermions at rest and has Sign[n] = 1. In 
the second configuration two fermions interchange their positions with one fermion 
stepping over the spatial boundary. If one uses periodic spatial boundary conditions 
this configuration has Sign[r2] = —1. Note that the same configuration would have 
Sign[n] = 1 when antiperiodic boundary conditions are used. 

The expectation value of a fermionic observable 0[n\ is given by 

(O)/ = ^E0NSign[n]exp(-5M). (2.14) 

Quantities of physical interest are the chiral condensate 

m[n] ='Y.{-lf^+^^+^^[n{x,t) - i), (2.15) 

and the corresponding chiral susceptibility 

X = ^{mf)f. (2.16) 

Up to now we have derived a path integral representation for the fermion system 
in terms of bosonic occupation numbers and a fermion sign factor that encodes Fermi 
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Figure 1: Two configurations of fermion occupation numbers in (1 + 1) dimensions. 
The shaded plaquettes carry the interaction. The dots represent occupied sites. In 
the second configuration two fermions interchange their positions. With periodic 
spatial boundary conditions this configuration has Sign[n] = — 1. 



statistics. The system without the sign factor is bosonic and is characterized by the 
positive Boltzmann factor exp(— S'[n]). Here the bosonic model is a quantum spin 
system with the Hamiltonian 



H = EiSlSl 



^3c3 



x+i 



Q^ C^ j_ r^ Co Co 

X x+i X x+l' 



(2.17) 



where S], = }^a\ is a spin 1/2 operator associated with the lattice site x that was 
labeled by /. The case G = 1 corresponds to the isotropic antiferromagnetic quan- 
tum Heisenberg model, G = represents the quantum XY- model, and G = —1 
corresponds to an isotropic ferromagnet. In the language of the spin model, the 
chiral condensate turns into the staggered magnetization 



mm 



E( 

x,t 



\Xi+X2+X;i q3 



(2.18) 



For sufficiently large G, for example in an antiferromagnet with G = 1, the staggered 
magnetization gets a non-zero expectation value at sufficiently low temperature, thus 
breaking the bosonic analog of chiral symmetry. It will turn out that the fermion 
sign does not change this behavior, and indeed chiral symmetry is spontaneously 
broken in the fermionic model as well. 



3 The Meron- Cluster Algorithm 



Let us first discuss the nature of the fermion sign problem. The fermionic path 
integral Zf = X^n Sign[n] exp(— 5'[n]) includes the fermion sign factor Sign[n] = ±1 



as well as a positive Boltzmann factor exp{—S[n]) that contains the action S[n] of 
the corresponding bosonic model with partition function Z^ = I]nexp(— S'[n]). A 
fermionic observable 0[n] is obtained in a simulation of the bosonic ensemble as 

iO)f = ^T.O[n]Sign[n]exp{-S[n]) = ^^^. (3.1) 

The average sign in the simulated bosonic ensemble is given by 

(Sign) = ^Yl SignW exp(-^[n]) = ^ = exp{-pVAf). (3.2) 

The last equality points to the heart of the sign problem. The expectation value of 
the sign is exponentially small in both the volume V and the inverse temperature 
P because the difference between the free energy densities A/ = ff — ft of the 
fermionic and bosonic systems is necessarily positive. 

Even in an ideal simulation of the bosonic ensemble which generates A^ com- 
pletely uncorrelated configurations, the relative statistical error of the sign is 



ASign _ V(Sign') - (Sign)^ _ expi/3VAf)_ 



(Sign) ViV(Sign) VN 

Here we have used Sign^ = 1. To determine the average sign with sufficient accuracy 
one needs to generate on the order of A^ = exp{2f3VAf) configurations. For large 
volumes and small temperatures this is impossible in practice. It is possible to 
solve one half of the problem if one can match any contribution —1 with another 
contribution 1 to give 0, such that only a few unmatched contributions 1 remain. 
Then effectively Sign = 0, 1 and hence Sign^ = Sign. This reduces the relative error 
to 



ASign ^ ^/(Sign) - (Sign)^ ^ exp{f3VAf/2) 

(Sign) ~ ^/N>{Sign) ~ y/W 

One gains an exponential factor in statistics, but one still needs to generate A^' = 
yN = exp{pVAf) independent configurations in order to accurately determine 
the average sign.Q This is because one generates exponentially many vanishing 
contributions before one encounters a contribution 1. As explained below, in our 
cluster algorithm an explicit matching of contributions —1 and 1 is achieved using 
an improved estimator. This solves one half of the sign problem. In a second step 
involving a Metropolis decision, our algorithm ensures that contributions and 1 
occur with similar probabilities. This saves another exponential factor in statistics 
and solves the other half of the sign problem. 

The meron-cluster fermion algorithm is based on a cluster algorithm for the cor- 
responding bosonic model without the sign factor. Bosonic quantum spin systems 

^The fact that an improved estimator alone cannot solve the sign problem was pointed out to 
one of the authors by H. G. Evertz a long time ago. 

8 



can be simulated very efficiently with cluster algorithms ||§, 0) §|- The ffist cluster 
algorithm for lattice fermions was described in 0. These algorithms can be imple- 



mented directly in the Euclidean time continuum ||T0[, i.e. the Suzuki- Trotter dis- 
cretization is not even necessary. The same is true for the meron-cluster algorithm. 
Here we discuss the algorithm for discrete time. The idea behind the algorithm is to 
decompose a configuration into clusters which can be flipped independently. Each 
lattice site belongs to exactly one cluster. When the cluster is flipped, the occupa- 
tion number of all the sites on the cluster is changed from n{x, t) to 1 — n{x, t), i.e. 
a cluster flip turns occupied into empty sites and vice versa. The decomposition of 
the lattice into clusters results from connecting neighboring sites on each individual 
space-time interaction plaquette following probabilistic cluster rules. A sequence of 
connected sites defines a cluster. In this case the clusters are sets of closed loops. 
The cluster rules are constructed to obey detailed balance. To show this we ffist 
write the plaquette Boltzmann factors as 

exp{—s[n{x, t),n{x + i,t), n{x, t + l),n{x + i,t + 1)]) = 

^^n{x,t),n{x,t+l)On(x+i,t),n{x+~i,t+l) + ^^n{x,t),l-n{x+t,t)^n{x,t+l),l-n{x+i,t+l) + 
^"n{x,t),n(x,t+'i-)"n(x+i,t),n{x+i,t+l)"n{x,t),l-n(x+i,t) "•" ^"n{x,t),n{x+i,t+l)"n{x+i,t),n{x,t+l) "•" 
ESn{x,t),n{x,t+l)Sn{x+lt),n{=^+i,t+l)^n{x,t),n{x+i,t)- (3-5) 

The various (5-functions specify which sites are connected and thus belong to the 
same cluster. The quantities A,B,...,E determine the relative probabilities for 
different cluster break-ups of an interaction plaquette. We only allow break-ups 
which generate legal configurations under cluster ffips. For example, A determines 
the probability with which sites are connected with their time-like neighbors, while 
B and D determine the probabilities for connections with space-like or diagonal 
neighbors, respectively. The quantities C and E determine the probabilities to put 
all four sites of a plaquette into the same cluster. This is possible for plaquette 
configurations [0,1,0,1] or [1,0,1,0] with a probability proportional to C and for 
configurations [0,0,0,0] or [1,1,1,1] with a probability proportional to E. The 
cluster rules are illustrated in table 1. 



Inserting the expressions from eq.( p.l2|) one finds 



eG 

exp(-s[0, 0, 0, 0]) = exp(-s[l, 1, 1,1]) = exp(-— ) = A + D + E, 

exp(-s[0, 1, 0, 1]) = exp(-s[l, 0, 1,0]) = cosh ^ = A + B + C, 
exp(-s[0, 1, 1,0]) = exp(-s[l, 0, 0, 1]) = sinh ^ = B + D. (3.6) 

For example, the probability to connect the sites with their time-like neighbors on 
a plaquette with configuration [0, 0, 0, 0] or [1, 1, 1, 1] is A/{A + D + E), while the 
probability for a connection with their diagonal neighbor is D/{A + D + E). All sites 
on such a plaquette are put into the same cluster with probability E/{A + D + E). 
Similarly, the probability for connecting space-like neighbors on a plaquette with 
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Table 1: Cluster break-ups of various plaquette configurations together with their 
relative probabilities A,B, ...,E. The dots represent occupied sites and the fat lines 
are the cluster connections. 



configuration [0,1,1,0] or [1,0,0,1] is B/{B 
connections is D/{B + D). 



D) and the probabihty for diagonal 



Eq. ( p. 51) can be viewed as a representation of the original model as a random clus- 
ter model. The cluster algorithm operates in two steps. First, a cluster break-up is 
chosen for each space-time interaction plaquette according to the above probabilities. 
This effectively replaces the original Boltzmann weight of a plaquette configuration 
with a set of constraints represented by the 5-functions associated with the chosen 
break-up. The constraints imply that occupation numbers of connected sites can 
only be changed together. In the second step of the algorithm every cluster is fiipped 
with probability 1/2. When a cluster is fiipped the occupation numbers of all sites 
that belong to the cluster are changed. Eq. p.(j| ) ensures that the cluster algorithm 
obeys detailed balance. To determine A,B,...,E we distinguish three cases. For 
G > 1 we solve eq.(|3.6|) by 
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(3.7) 



G = 0, 



(3. 
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and, finally, for G < — 1 

A = cosh |, B = C = 0, D = sinh L E = exp(-^) - exp(|). (3.9) 

As an example, let us consider the antiferromagnetic quantum Heisenberg model, 
i.e. G = 1, and hence 

A = exp(--), B = sinh i, C = D = E = 0. (3.10) 

Consequently, on plaquette configurations [0, 0, 0, 0] or [1, 1, 1, 1] one always chooses 
time-like connections between sites, and for configurations [0, 1, 1, 0] or [1, 0, 0, 1] one 
always chooses space-like connections. For configurations [0, 1, 0, 1] or [1, 0, 1, 0] one 
chooses time-like connections with probability p = A/ [A + B) = 2/[l + exp(e/2)] 
and space-like connections with probability 1 — p = B/{A + B). Indeed, this is the 
algorithm that was used in 0]. It is extremely efficient, has almost no detectable 
autocorrelations, and its dynamical exponent for critical slowing down is compatible 
with zero. 

Let us now consider the effect of a cluster flip on the fermion sign. It is obvious 
that the flip of a cluster either leads to a sign change or it leaves the sign unchanged. 
In general, the effect of the flip of a specific cluster on the fermion sign depends on 
the orientation of the other clusters. For example, a cluster whose flip does not 
change the sign now, may very well change the sign after other clusters have been 
flipped. In other words, the clusters affect each other in their effect on the fermion 
sign. This makes it very difficult to understand the effect of the various cluster flips 
on the topology of the fermion world-lines and thus on Sign[n]. As a consequence, 
for the most general model described above, we don't know how to solve the fermion 
sign problem. However, there are clusters whose effect on Sign[n] is independent of 
the orientation of all the other clusters. Specifically, these are those clusters that 
do not contain any diagonal break-ups. To ensure that no clusters contain such 
break-ups we limit ourselves to models which have D = 0. According to eq.( p.7|) 
this is the case for G > 1. 

Once we forbid diagonal cluster break-ups, i.e. when D = 0, the clusters have 
a remarkable property with far reaching consequences: each cluster can then be 
characterized by its effect on the fermion sign independent of the orientation of all 
the other clusters. We refer to clusters whose flip changes Sign[n] as merons, while 
clusters whose flip leaves Sign[n] unchanged are called non- merons. The flip of a 
meron-cluster permutes the fermions and changes the topology of the fermion world- 
lines. The term meron has been used before to denote half-instantons. For example, 
the meron-clusters in the algorithm for the 2-d 0(3) model at non- zero vacuum 



angle 9 are indeed half-instantons ||I2[. A configuration with an odd permutation of 
fermions has Sign[r2] = — 1 and can be viewed as an instanton, while a configuration 
with an even permutation of fermions is topologically trivial and has Sign[ra] = 
1. The flip of a meron-cluster changes an instanton into a topologically trivial 
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configuration and therefore is a lialf-instanton. In (2 + 1) dimensions particles can 
liave any statistics. Tlie anyon statistics is cliaracterized by a pliase factor expliOH), 
where H is the integer- valued Hopf number. In that case, the flip of a meron-cluster 
changes the Hopf number by one, and the meron is indeed a half-Hopf-instanton. 
The number of merons in a configuration is always even. This follows when one 
considers flipping all clusters, and thus changing the occupation of all lattice sites, 
which leaves the fermion sign unchanged. If the number of merons would be odd, 
flipping all clusters would change the fermion sign, which implies that the number 
of merons must be even. An example of a meron-cluster is given in figure 2, which 
contains the same fermion configurations as figure 1. When the meron-cluster is 
flipped the first configuration with Sign[n] = 1 turns into the second configuration 
with Sign[n] = — 1. 
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Figure 2: The same fermion configurations as in figure 1 together with a meron- 
cluster represented by the fat line. The other clusters are not shown. Flipping the 
meron-cluster changes one configuration into the other and changes the fermion 
sign. 

The meron concept alone allows us to gain an exponential factor in statistics. 
Since all clusters can be flipped independently, one can construct an improved es- 
timator for Sign[n] by averaging analytically over the 2^^ configurations obtained 
by flipping the Nc clusters in the configuration in all possible ways. For config- 
urations that contain merons the average Sign[n] is zero because flipping a single 
meron leads to a cancellation of contributions ±1. Hence, only the configurations 
without merons contribute to Sign[n]. The vast majority of configurations contains 
merons and now contributes an exact to Sign[n] instead of a statistical average of 
contributions ±1. In this way the improved estimator leads to an exponential gain 
in statistics. 

Still, as it stands it is not guaranteed that the contributions from the zero-meron 
sector will always be positive. With no merons in the configuration it is clear that 
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the fermion sign remains unchanged under cluster flip, but one could still have 
Sign[n] = — 1. Fortunately, there is a way to guarantee that Sign[n] = 1 in the zero- 
meron sector. In that case the contributions to Sign[n] are from configurations 
containing meron-clusters and 1 from configurations without merons. According to 
the previous discussion, this solves one half of the fermion sign problem. Let us 
describe how one can make sure that a configuration without merons always has 
Sign[?T,] = 1. This is possible when one does not allow the cluster break-ups charac- 
terized by the amplitudes D and E in eq.(|3.5|). This again limits us to models with 



G > 1, for which we have D = E = according to eq.( p.7|) . The remaining cluster 
break-ups with amplitudes A, B, C have a very important property. They guarantee 
that sites inside a cluster obey a pattern of staggered occupation, i.e. either the 
even sites (with xi -|- X2 + x^ even) within the cluster are all occupied and the odd 
sites are all empty, or the even sites are all empty and the odd sites are all occupied. 
This guarantees that the clusters can be flipped such that one reaches the totally 
staggered reference configuration in which at all times all even sites are occupied 
and all odd sites are empty. In the corresponding antiferromagnet this configura- 
tion represents a completely ordered state with a staggered magnetization. In the 
half-occupied reference configuration (which is the first configuration in figure 1) all 
fermions are at rest, no fermions are permuted during the Euclidean time evolution, 
and thus Sign[n] = 1. Since any configuration can be turned into the reference 
configuration by appropriate cluster flips, this is particularly true for configurations 
without merons. Since the totally staggered configuration has Sign[n] = 1 and the 
fermion sign remains unchanged when a non-meron-cluster is flipped, all configura- 
tions without merons have Sign[n] = 1. Instead of a sequence of ±1 for Sign[n], we 
now have contributions and 1. As discussed before, this only solves one half of the 
fermion sign problem. Before we can solve the other half of the problem we must 
discuss improved estimators for the physical observables. 

Let us consider an improved estimator for (\l/\l/[?2])^Sign[n] which is needed to de- 
termine the chiral susceptibility x- The total chiral condensate, \l/\l/[n] = "^Zc^^Ci 
is a sum of cluster contributions 

^^c = ^ E {-ir^^'^''Kn{x^)-\). (3.11) 

" {x,t)&C ^ 

When a cluster is flipped, its condensate contribution changes sign. In a config- 
uration without merons, where Sign[?7.] = 1 for all relative cluster flips, the aver- 
age of (\l/^[n])^Sign[?T,] over all 2^"^ configurations is Z^cl^^cP- For configura- 
tions with two merons the average is 2|\I'\I^cJ|\I'\I'c72l where C\ and C2 are the two 
meron-clusters. Configurations with more than two merons do not contribute to 
(^\l/[?2])^Sign[n]. The improved estimator for the susceptibility is hence given by 

where N is the number of meron-clusters in a configuration. Thus, to determine x 
one must only sample the zero- and two-meron sectors. 
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The probability to find a configuration without merons is exponentially small in 
the space-time volume since it is equal to (Sign). Thus, although we have increased 
the statistics tremendously with the improved estimators, without a second step 
one would still need an exponentially large statistics to accurately determine x- 
Fortunately, the numerator in equation (|3.12|) receives contributions from the zero- 
and two-meron sectors only, while the denominator gets contributions only from 
the zero-meron sector. One can hence restrict oneself to the zero- and two-meron 
sectors and never generate configurations with more than two merons. This enhances 
both the numerator and the denominator by a factor that is exponentially large in 
the volume, but leaves the ratio of the two invariant. One purpose of the second 
step of the meron-cluster algorithm is to eliminate all configurations with more 
than two merons. To achieve this, we start with an initial configuration with zero 
or two merons. For example, a completely occupied configuration has no merons. 
We then visit all plaquette interactions one after the other and choose new pair 
connections between the four sites according to the above cluster rules. If the new 
connection increases the number of merons beyond two, it is not accepted and the old 
connection is kept for that plaquette. This procedure obeys detailed balance because 
configurations with more than two merons do not contribute to the observable we 
consider. This simple reject step eliminates almost all configurations with weight 
and is the essential step to solve the other half of the fermion sign problem. 

Assuming a dilute gas of meron and non-meron-clusters of typical space-time 
volume \C\ one finds a ratio p(0)/p(2) oc {\C\/V P)"^ of the probabilities p(0) andp(2) 
to have zero or two merons. Hence, as long as the cluster size does not grow with the 
space-time volume, most configurations would have two merons and therefore would 
still have weight 0. Without further improvements, one would still need statistics 
quadratic (but no longer exponential) in the space-time volume to get an accurate 
average sign. The remaining problem can be solved with a reweighting technique 
similar to the one used in [|I^. To enhance the zero-meron configurations in a 
controlled way, we introduce trial probabilities pt(0) and Pt(2) which determine the 
relative weight of the zero- and two-meron sector. The trial distribution pt{N) for 
A^ > 2 is set to infinity. The distribution pt{N) is used in a Metropolis accept-reject 
step for the newly proposed pair connection on a specific plaquette interaction. A 
new pair connection that changes the meron number from A^ to A^' is accepted 
with probability p = inm[l, pt{N) / pt{N')]. In particular, configurations with A^' > 
2 are never generated because then pt{N') = cxd and p = 0. After visiting all 
plaquette interactions, each cluster is flipped with probability 1/2 which completes 
one update sweep. After reweighting, the zero- and two-meron configurations appear 
with similar probabilities. This completes the second step in our solution of the 
fermion sign problem. The reweighting of the zero- and two-meron configurations is 
taken into account in the final expression for the chiral susceptibility as 

^ (Ec |^^c|''^JV,o PtjO) + 2\^^cA\^^c,\Sn,2 Pt{2)) 
^ VP{5r,,o pM) • ^ ■ ' 
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L 


/3 


(Sign) 


vm/vm 


(Sign)^ 


X 


4 


0.6 


0.838(1) 


0.5/0.5 


0.845(1) 


0.554(1) 


4 


0.7 


0.710(3) 


0.5/0.5 


0.726(2) 


0.936(2) 


4 


0.8 


0.534(4) 


0.5/0.5 


0.566(2) 


1.678(5) 


4 


0.9 


0.357(3) 


0.5/0.5 


0.405(2) 


3.13(1) 


4 


0.948 


0.282(2) 


0.3/0.7 


0.537(3) 


4.22(2) 


6 


0.948 


0.0556(7) 


0.2/0.8 


0.398(3) 


9.83(8) 


8 


0.948 


0.0020(4) 


0.1/0.9 


0.361(8) 


16.6(5) 


10 


0.948 


— 


0.1/0.9 


0.178(3) 


26.5(8) 


12 


0.948 


— 


0.05/0.95 


0.17(1) 


37(1) 


14 


0.948 


— 


0.02/0.98 


0.20(1) 


51(1) 


16 


0.948 


— 


0.01/0.99 


0.22(2) 


65(2) 



Table 2: Numerical results for the non-reweighted (Sign), the reweighted {Sign)r 
and X obtained with a reweighting factor pt{0)/pt{2) on lattices of spatial size L at 
inverse temperature f3. For the larger volumes the non-reweighted (Sign) is too small 
to be measured. 

The optimal ratio pt{0)/pt{2), which minimizes the statistical error of X; can be 
estimated by gradually increasing the volume and the inverse temperature to their 
desired values. 



4 Numerical Results 



We have simulated the staggered fermion model with G = 1 on antiperiodic spatial 
volumes L^ with L = 4,6, ..., 16 and at various inverse temperatures (3 G [0.5, 1.2] 
which includes the critical point. In the Euclidean time direction we have used 
M = 4, i.e. 24 time-slices. In all cases, we have performed at least 1000 thermaliza- 
tion sweeps followed by 10000 measurements. One sweep consists of a new choice 
of the cluster connections on each interaction plaquette and a flip of each cluster 
with probability 1/2. To estimate an optimal ratio pt{0)/pt{2) of the reweighting 
probabilities, we have first performed runs without reweighting. The observed rel- 
ative probabilities of the zero- and two-meron sectors were used as an estimate for 
Pt{0)/pt{2) in production runs. Especially for the larger volumes, reweighting is 
necessary in order to obtain accurate results. 

Some of our data for the chiral susceptibility x ^-^e contained in table 2. The 
table also includes the non-reweighted (Sign) , the value of the used reweighting fac- 
tor pt{0)/pt{2) as well as the reweighted (Sign),,. Note that (Sign)^ is the fraction 
of zero-meron configurations that the algorithm generates by sampling the zero- 
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and two-meron sectors only. This quantity is typically a lot bigger than the orig- 
inal non-reweighted (Sign), which is the fraction of zero-meron configurations in 
the space of all configurations including those with many merons. In particular, 
in large space-time volumes the staggered fermion model suffers from a very severe 
sign problem. Figure 3 shows the probability to have a certain number of merons 
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Figure 3: The probability of having a certain number of merons for various spatial 
sizes L = 4, 6, 8 and 10 at p = 0.948. 



in an algorithm that samples all meron-sectors without reweighting. For small vol- 
umes the zero-meron sector and hence (Sign) are relatively large, while multi-meron 
configurations are rare. On the other hand, in larger volumes the vast majority of 
configurations has a large number of merons and hence (Sign) is exponentially small. 
For example, an extrapolation from smaller volumes gives a rough estimate for the 
non-reweighted (Sign) ^ 10^^° on the 16^ lattice aX (3 = 0.948, while the reweighted 
(Sign)r = 0.22(2). Hence, to achieve a similar accuracy without the meron-cluster 
algorithm one would have to increase the statistics by a factor 10^°, which is obvi- 
ously impossible in practice. In fact, at present there is no other method that can 
be used to simulate this model. 

Figure 4 shows the chiral susceptibility x ^-s a function of f3 for various spatial 
sizes L. At high temperatures (small f3) x is almost independent of the volume, in- 
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dicating that chiral symmetry is unbroken. On the other hand, at low temperatures 
X increases with the volume, which implies that chiral symmetry is spontaneously 
broken. To study the critical behavior in detail, we have performed a finite-size 
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Figure 4: The chiral susceptibility x o,s a function of the inverse temperature (3 for 
various spatial sizes L = 4, 6, 8 and 10. 



scaling analysis for x focusing on a narrow range (3 G [0.9, 0.98] around the critical 
point. Since a Z(2) chiral symmetry gets spontaneously broken at finite temperature 
in this (3 + l)-d model, one expects to find the critical behavior of the 3-d Ising 



model. The corresponding finite-size scaling formula valid close to /3c is [0 



x(L,/5)=a(x)+%)L^^ 

a{x) = ao + flix + a2x'^ + ..., x = P — Pc, 

b{y) = bo + hy + b,y' + ...^y={(3- (3^)L^/\ 



(4.1) 



For the 3-d Ising model the critical exponents are given by z/ = 0.630(1) and 7/// = 
1.963(3) m. Fitting our data, we find u = 0.63(4) and 7/// = 1.98(2), which 
indicates that the chiral transition of the staggered fermion model is indeed in the 
3-d Ising universality class. The fit gives (3c = 0.948(3). In figure 5 we have taken 
the values of the critical exponents from the 3-d Ising model, and we have plotted 
XJU^I" as a function oi y = {/3 — /3c)L^^'^- For large enough L one can neglect the 
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term a(x) in eq.( |4.1| ) and one obtains x/L^^'^ = Kv)- We have varied the value 
of f3c and found that indeed all data can be collapsed on one universal curve. The 
resulting critical inverse temperature is (3c = 0.948(3) in agreement with the previous 
fit. At the estimated value of /3c, we have performed simulations on larger spatial 
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Figure 5: Finite-size scaling behavior of the chiral susceptibility x- The data for 
various spatial sizes L = 4, 6, 8 and 10 fall on one universal curve. 

volumes up to 16^. The results for x are shown in figure 6 together with a fit that 
gives an independent estimate of 7/z/ = 1.98(2). All of this supports the claim that 
the chiral transition in our model is Ising-like. 



5 Conclusions 



We have applied a new fermion simulation technique — the meron-cluster algorithm 
— to a model of staggered fermions. In contrast to standard methods which inte- 
grate out the fermions and are left with a non-local bosonic action, we describe the 
fermions in a Fock state occupation number basis and thus keep a local bosonic 
action. The fermion permutation sign arises due to the Pauli exclusion principle as 
a non-local factor associated with non-trivial topology of the fermion world-lines. 
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Figure 6: Finite-size scaling behavior of the chiral susceptibility x (is a function of 
the spatial size L at the estimated critical inverse temperature Pc = 0.948. A fit of 
the volume-dependence (solid line) gives the critical exponent 7/z/ = 1.98(2) of the 
3-d Ising model. 



The sign leads to severe cancellations which makes standard simulation techniques 
impossible to use. The decomposition of a configuration into clusters allows us 
to disentangle the complicated topology of the fermion world-lines. In particular, 
a meron-cluster identifies a pair of configurations with equal weight and opposite 
sign. This results in an explicit cancellation of two contributions ±1 to the path 
integral, such that only the zero-meron sector contributes to the partition function. 
Observables like the chiral susceptibility x receive contributions only from the zero- 
and two-meron sectors. To measure x oiie can hence eliminate all sectors with more 
than two merons, which leads to an exponential gain in statistics and to a complete 
solution of the fermion sign problem. 

The meron-cluster algorithm allowed us to simulate a staggered fermion model 
for which standard fermion methods suffer from a severe sign problem. The model 
has two fiavors and a Z(2) chiral symmetry which is spontaneously broken at low 
temperatures. Applying finite-size scaling methods to the high-precision numerical 
data for the chiral susceptibility, we extracted critical exponents compatible with 
those of the 3-d Ising model. This is the expected behavior based on universality 
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arguments and dimensional reduction. It would be interesting to apply our method 
to (2+1) dimensions. The A^ = 4 flavor case was studied in p| and it was verified that 
the model is in the 2-d Ising universality class. The standard fermion algorithm that 
was used in that study does not work for A^ < 4 due to the fermion sign problem. 
The meron-cluster algorithm solves the sign problem and can be used to explore 
those models. 

Even in cases without a sign problem, the meron-cluster algorithm is more ef- 
ficient than standard fermion simulation methods. However, one should keep in 
mind that the meron-cluster algorithm is not always applicable. For example, the 
meron concept applies only when the clusters are independent in their effect on the 
fermion sign. In addition, it must always be possible to fiip the clusters such that 
one reaches a reference configuration with a positive sign. Otherwise, some contri- 
butions from the zero-meron sector could be negative. For example, in our model 
these restrictions led to G > 1, i.e. to a sufficiently strong four- fermion coupling. In 
this paper, we have not attempted to take the continuum limit of the lattice theory. 
Instead, we have studied the chiral phase transition at a finite lattice spacing corre- 
sponding to G = 1. The restriction G > 1 would prevent us from approaching the 
continuum limit if it corresponds to Gc < 1. It should be pointed out that we can 
study the universal behavior of the chiral transition without taking the continuum 
limit. Of course, it is not excluded that appropriate modifications of the algorithm 
might allow us to work at G < 1. 

A natural arena for applications of the meron-cluster algorithm are quantum 
link models |T^ which are used in the D-theory formulation of QCD |1^, |16|, [1^ . In 



D-theory continuous classical fields arise from the dimensional reduction of discrete 
quantum variables. In these models the continuum limit is approached by vary- 
ing the extent of an extra dimension, not by varying the value of a bare coupling 
constant. Hence, restrictions such as G > 1 would not prevent us from taking the 
continuum limit. In quantum link QCD the quarks arise as domain wall fermions. 
The application of meron-cluster algorithms to domain wall fermions is in progress. 
Also there are many applications to sign problems in condensed matter physics. In- 
vestigations of antiferromagnets in a magnetic field and of systems in the Hubbard 
model family will be reported elsewhere. 

At present, the meron-cluster algorithm is the only method that allows us to solve 
the fermion sign problem. A severe sign problem arises in lattice QCD calculations 
at non-zero baryon number due to a complex action. It is therefore natural to ask 
if our algorithm can be applied to this case. At non-zero chemical potential the 2-d 
0(3) model, which is a toy model for QCD, also suffers from a sign problem due 
to a complex action. When applied to the D-theory formulation of this model, the 
meron-cluster algorithm solves the sign problem completely. It remains to be seen 
if a similar result can be achieved for QCD. 
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